## Read data	
Data <- read.table("HMS_Data.raw", header=T)

## Compute observed strata probabilities

# Question 3
obs_Q3 <- ifelse(is.na(Data$CQ3a_Clinic_or_Hospitals)==F & is.na(Data$FQ3a_Clinic)==F, 1, 0)
Q3_total <- subset(Data, obs_Q3 == 1)
Q3_00 <- ifelse(Q3_total$CQ3a_Clinic_or_Hospitals == 0 & Q3_total$FQ3a_Clinic ==0, 1, 0)
Q3_01 <- ifelse(Q3_total$CQ3a_Clinic_or_Hospitals == 0 & Q3_total$FQ3a_Clinic ==1, 1, 0)
Q3_10 <- ifelse(Q3_total$CQ3a_Clinic_or_Hospitals == 1 & Q3_total$FQ3a_Clinic ==0, 1, 0)Q3_11 <- ifelse(Q3_total$CQ3a_Clinic_or_Hospitals == 1 & Q3_total$FQ3a_Clinic ==1, 1, 0)
Pr00_Q3 <- sum(Q3_00)/sum(obs_Q3)
Pr01_Q3 <- sum(Q3_01)/sum(obs_Q3)
Pr10_Q3 <- sum(Q3_10)/sum(obs_Q3)
Pr11_Q3 <- sum(Q3_11)/sum(obs_Q3)
P.3 <- c(Pr00_Q3, Pr01_Q3, Pr10_Q3, Pr11_Q3)


# Question 4a
obs_Q4a <- ifelse(is.na(Data$CQ4a_Ed_Prim_or_Secondary)==F & is.na(Data$FQ4a_Ed_Prim)==F, 1, 0)
Q4a_total <- subset(Data, obs_Q4a == 1)
Q4a_00 <- ifelse(Q4a_total$CQ4a_Ed_Prim_or_Secondary == 0 & Q4a_total$FQ4a_Ed_Prim ==0, 1, 0)
Q4a_01 <- ifelse(Q4a_total$CQ4a_Ed_Prim_or_Secondary == 0 & Q4a_total$FQ4a_Ed_Prim ==1, 1, 0)
Q4a_10 <- ifelse(Q4a_total$CQ4a_Ed_Prim_or_Secondary == 1 & Q4a_total$FQ4a_Ed_Prim ==0, 1, 0)
Q4a_11 <- ifelse(Q4a_total$CQ4a_Ed_Prim_or_Secondary == 1 & Q4a_total$FQ4a_Ed_Prim ==1, 1, 0)
Pr00_Q4a <- sum(Q4a_00)/sum(obs_Q4a)
Pr01_Q4a <- sum(Q4a_01)/sum(obs_Q4a)
Pr10_Q4a <- sum(Q4a_10)/sum(obs_Q4a)
Pr11_Q4a <- sum(Q4a_11)/sum(obs_Q4a)
P.4a <- c(Pr00_Q4a, Pr01_Q4a, Pr10_Q4a, Pr11_Q4a)


# Question 7b
obs_Q7b <- ifelse(is.na(Data$CQ7b_Tr_Roads_or_Services)==F & is.na(Data$FQ7b_Tr_Roads)==F, 1, 0)
Q7b_total <- subset(Data, obs_Q7b == 1)
Q7b_00 <- ifelse(Q7b_total$CQ7b_Tr_Roads_or_Services == 0 & Q7b_total$FQ7b_Tr_Roads ==0, 1, 0)
Q7b_01 <- ifelse(Q7b_total$CQ7b_Tr_Roads_or_Services == 0 & Q7b_total$FQ7b_Tr_Roads ==1, 1, 0)
Q7b_10 <- ifelse(Q7b_total$CQ7b_Tr_Roads_or_Services == 1 & Q7b_total$FQ7b_Tr_Roads ==0, 1, 0)
Q7b_11 <- ifelse(Q7b_total$CQ7b_Tr_Roads_or_Services == 1 & Q7b_total$FQ7b_Tr_Roads ==1, 1, 0)
Pr00_Q7b <- sum(Q7b_00)/sum(obs_Q7b)
Pr01_Q7b <- sum(Q7b_01)/sum(obs_Q7b)
Pr10_Q7b <- sum(Q7b_10)/sum(obs_Q7b)
Pr11_Q7b <- sum(Q7b_11)/sum(obs_Q7b)
P.7b <- c(Pr00_Q7b, Pr01_Q7b, Pr10_Q7b, Pr11_Q7b)

# Question 11a
obs_Q11a <- ifelse(is.na(Data$CQ11a_Discount_Binary)==F & is.na(Data$FQ11a_Discount_Bin)==F, 1, 0)
Q11a_total <- subset(Data, obs_Q11a == 1)
Q11a_00 <- ifelse(Q11a_total$CQ11a_Discount_Binary == 0 & Q11a_total$FQ11a_Discount_Bin ==0, 1, 0)
Q11a_01 <- ifelse(Q11a_total$CQ11a_Discount_Binary == 0 & Q11a_total$FQ11a_Discount_Bin ==1, 1, 0)
Q11a_10 <- ifelse(Q11a_total$CQ11a_Discount_Binary == 1 & Q11a_total$FQ11a_Discount_Bin ==0, 1, 0)
Q11a_11 <- ifelse(Q11a_total$CQ11a_Discount_Binary == 1 & Q11a_total$FQ11a_Discount_Bin ==1, 1, 0)
Pr00_Q11a <- sum(Q11a_00)/sum(obs_Q11a)
Pr01_Q11a <- sum(Q11a_01)/sum(obs_Q11a)
Pr10_Q11a <- sum(Q11a_10)/sum(obs_Q11a)
Pr11_Q11a <- sum(Q11a_11)/sum(obs_Q11a)
P.11a <- c(Pr00_Q11a, Pr01_Q11a, Pr10_Q11a, Pr11_Q11a)


## Obtain upper and lower bounds
## read libraries
library(lpSolve)

## define functions 

sens.a13 <- function(P,dir,e,d){
evpoints <- round(1/e) + 1
corpoints <- 1/d
q <- 0
rho <- 0
bound <- rep(NA, evpoints)
pi.cc <- rep(NA, evpoints)
pi.ca <- rep(NA, evpoints)
pi.cn <- rep(NA, evpoints)
pi.cd <- rep(NA, evpoints)
pi.ac <- rep(NA, evpoints)
pi.aa <- rep(NA, evpoints)
pi.an <- rep(NA, evpoints)
pi.ad <- rep(NA, evpoints)
pi.nc <- rep(NA, evpoints)
pi.na <- rep(NA, evpoints)
pi.nn <- rep(NA, evpoints)
pi.nd <- rep(NA, evpoints)
pi.dc <- rep(NA, evpoints)
pi.da <- rep(NA, evpoints)
pi.dn <- rep(NA, evpoints)
pi.dd <- rep(NA, evpoints)
status <- rep(NA, evpoints)
Q <- rep(NA, evpoints)
Px0 <- P[1] + P[3]
Px1 <- P[2] + P[4]
sharpb <- rep(NA, corpoints)
corr <- rep(NA, corpoints)
for(j in 1:corpoints){
    rho <- (j-1)*d
    for(i in 1:evpoints){
        q <- (i-1)*e
        f.obj <- c(1,1,1,1,0,0,0,0,0,0,0,0,-1,-1,-1,-1)
        f.con <- matrix(c(1-q,0,1-q,0,0,0,0,0,1-q,0,1,q,0,0,q,q,
                0,1-q,0,1-q,0,0,0,0,q,1,0,1-q,q,q,0,0,
                0,0,q,q,1-q,0,1,q,0,0,0,0,1-q,0,1-q,0,
                q,q,0,0,q,1,0,1-q,0,0,0,0,0,1-q,0,1-q,
                0,1,1,2,0,1,1,2,0,1,1,2,0,1,1,2,
                1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
                0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
                0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,
                0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,
                0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,
                0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,
                0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,
                0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,
                0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,
                0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
                0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,
                0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,
                0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,
                0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,
                0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,
                0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1), nrow=21, byrow=T)
        f.dir <- c("=","=","=","=","<",
            ">=",">=",">=",">=",">=",">=",">=",">=",
            ">=",">=",">=",">=",">=",">=",">=",">=")
        f.rhs <- c(P,rho,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0)
        r <- lp(dir, f.obj, f.con, f.dir, f.rhs)
        bound[i] <- ifelse(r$status==2, NA, r$objval)
        pi.cc[i] <- r$solution[1]
        pi.ca[i] <- r$solution[2]
        pi.cn[i] <- r$solution[3]
        pi.cd[i] <- r$solution[4]
        pi.ac[i] <- r$solution[5]
        pi.aa[i] <- r$solution[6]
        pi.an[i] <- r$solution[7]
        pi.ad[i] <- r$solution[8]
        pi.nc[i] <- r$solution[9]
        pi.na[i] <- r$solution[10]
        pi.nn[i] <- r$solution[11]
        pi.nd[i] <- r$solution[12]
        pi.dc[i] <- r$solution[13]
        pi.da[i] <- r$solution[14]
        pi.dn[i] <- r$solution[15]
        pi.dd[i] <- r$solution[16]
        status[i] <- r$status
        Q[i] <- q
        }
    sharpb[j] <- ifelse(dir=="max", max(na.exclude(bound)), min(na.exclude(bound)))
    corr[j] <- rho
    }
return(cbind(corr, sharpb))
}

Qby <- 0.001
rhoby <- 0.001

sensresQ3.13.u <- sens.a13(P.3, "max", Qby, rhoby)
sensresQ3.13.l <- sens.a13(P.3, "min", Qby, rhoby)

sensresQ4a.13.u <- sens.a13(P.4a, "max", Qby, rhoby)
sensresQ4a.13.l <- sens.a13(P.4a, "min", Qby, rhoby)

sensresQ7b.13.u <- sens.a13(P.7b, "max", Qby, rhoby)
sensresQ7b.13.l <- sens.a13(P.7b, "min", Qby, rhoby)

sensresQ11a.13.u <- sens.a13(P.11a, "max", Qby, rhoby)
sensresQ11a.13.l <- sens.a13(P.11a, "min", Qby, rhoby)

# Create plots

pdf("fig3.pdf", width = 11, height = 8.5)
par(mfrow=c(2,2))

# Question 3
plot(c(0,1), c(-1,1), type="n", asp=0.3,
    main = "Q3: Clinics (0) or Hospitals (1)?",
    xlab = expression(rho),
    ylab = "Average Treatment Effect")
lines(sensresQ3.13.u[,1], sensresQ3.13.u[,2], col="darkblue", lty="solid", lwd=3)
lines(sensresQ3.13.l[,1], sensresQ3.13.l[,2], col="darkblue", lty="solid", lwd=3)
lines(c(0.142, 0.142), c(0, -1.1), lty="dashed", lwd=2)
abline(h = P.3[4]/(P.3[2]+P.3[4]) - P.3[3]/(P.3[1]+P.3[3]), lty="dotted", lwd=3)
abline(0,0)

# Question 4a
plot(c(0,1), c(-1,1), type="n", asp=0.3,
    main = "Q4a: Primary (0) or Secondary (1) Education?",
    xlab = expression(rho),
    ylab = "Average Treatment Effect")
lines(sensresQ4a.13.u[,1], sensresQ4a.13.u[,2], col="darkblue", lty="solid", lwd=3)
lines(sensresQ4a.13.l[,1], sensresQ4a.13.l[,2], col="darkblue", lty="solid", lwd=3)
lines(c(0.246, 0.246), c(0, -1.1), lty="dashed", lwd=2)
abline(h = P.4a[4]/(P.4a[2]+P.4a[4]) - P.4a[3]/(P.4a[1]+P.4a[3]), lty="dotted", lwd=3)
abline(0,0)

# Question 7b
plot(c(0,1), c(-1,1), type="n", asp=0.3,
    main = "Q7b: Roads (0) or Public Transportation (1)?",
    xlab = expression(rho),
    ylab = "Average Treatment Effect")
lines(sensresQ7b.13.u[,1], sensresQ7b.13.u[,2], col="darkblue", lty="solid", lwd=3)
lines(sensresQ7b.13.l[,1], sensresQ7b.13.l[,2], col="darkblue", lty="solid", lwd=3)
lines(c(0.005, 0.005), c(0, -1.1), lty="dashed", lwd=2)
abline(h = P.7b[4]/(P.7b[2]+P.7b[4]) - P.7b[3]/(P.7b[1]+P.7b[3]), lty="dotted", lwd=3)
abline(0,0)

# Question 11a
plot(c(0,1), c(-1,1), type="n", asp=0.3,
    main = "Q11a: Consume (0) or Invest (1) Windfalls?",
    xlab = expression(rho),
    ylab = "Average Treatment Effect")
lines(sensresQ11a.13.u[,1], sensresQ11a.13.u[,2], col="darkblue", lty="solid", lwd=3)
lines(sensresQ11a.13.l[,1], sensresQ11a.13.l[,2], col="darkblue", lty="solid", lwd=3)
lines(c(0.056, 0.056), c(0, -1.1), lty="dashed", lwd=2)
abline(h = P.11a[4]/(P.11a[2]+P.11a[4]) - P.11a[3]/(P.11a[1]+P.11a[3]), lty="dotted", lwd=3)
abline(0,0)

dev.off()

